Speckle instability: coherent effects in nonlinear disordered media 
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We numerically investigate the properties of speckle patterns formed by nonlinear point scatterers. 
We show that, in the weak localization regime, dynamical instability appears, eventually leading to 
chaotic behavior of the system. Analysing the statistical properties of the instability thresholds for 
different values of the system size and disorder strength, a scaling law is emphasized. The later is 
also found to govern the smallest decay rate of the linear system, putting thus forward the crucial 
importance of interference effects. This is also underlined by the fact that coherent backscattering 
is still observed even in the chaotic regime. 
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As first described by Anderson in 1958, the impact 
of disorder on the transport properties of waves, de- 
pending on the dimensionality and the disorder strength, 
ranges from weak to strong localization. In the case 
of matter waves, the localization of Bose- Einstein con- 
densates (BEG) is, at present, a very active research 
topic investigated by several experimental and theoretical 
groups [Tl[51[31lll[51[Sl[71[Sj, in particular with the recent 
observation of the (ID) localization of matter waves in a 
disordered optical potentials [9l[10]. Although, in these 
experiments, the atom-atom interactions were negligible, 
the basic question remains how effects of interference be- 
tween multiply scattered waves, such as weak or strong 
localization, are affected by interactions. 

BEG's in disordered potentials appear to be good can- 
didates to study these questions, since, in the mean field 
regime, the condensate is still described by a single co- 
herent wave function governed by a nonlinear wave equa- 
tion (Gross-Pitaevski equation). Thus, the condensate in 
principle retains its ability to display interference effects 
also in presence of (not too strong) interactions. Similar 
nonlinear equations describe propagation of light in dis- 
ordered nonlinear media In contrast, the situation 
is quite different in the case of electronic transport [T2], 
where the interactions combined with finite temperature 
effects give rise to dephasing, which in general destroy 
the disorder-induced coherent effects. 

Even if the theoretical description of coherent effects in 
nonlinear disordered systems is far from complete, an im- 
portant step was done recently by the development of a 
diagrammatic theory for coherent backscattering in pres- 
ence of nonlinearity [13]. This approach relies on the as- 
sumption of a unique stationary solution of the nonlinear 
wave equation under consideration. This assumption is 
expected to be valid for very weak nonlinearities (which, 
however, may still considerably affect the height of the 
coherent backscattering cone [HI US])- On the other hand. 



it is known that larger nonlinearities can induce speckle 
instabilities, such that no stationary state is reached at 
long times [SJ [T3] . A clear explanation for the physical 
origin of this effect, however, is still missing. Theoretical 
predictions of the nonlinearity threshold above which in- 
stabilities develop, were attempted in |1HI15]. but, to our 
knowledge, these results have not yet been confirmed by 
experiments or numerical simulations. For this reason, 
we investigate, in this letter, the dynamical properties of 
the speckle patterns generated by nonlinear point scat- 
terers. 

More precisely, we consider an assembly of N point- 
like scatterers located at randomly chosen positions r^, 
i — I, . . . , N, inside a sphere of volume V illuminated 
by a plane wave k l ■ The field amplitude Ei at scatterer 
i results as the sum of the incident wave and the field 
radiated by all other scatterers: 



E 



(1) 



where k = Ik^l = 27r/A, and the field amplitudes are 
measured in units of the incident plane wave amplitude. 
dj describes the dipole induced inside scatterer j. For 
simplicity, we will consider only scalar fields in this letter. 

As for the induced dipole, we consider a nonlinear re- 
sponse d = g (l-Ep) E, with 



(2) 



For a = 0, this model corresponds to a resonant point 
scatterer. The nonlinearity a induces a phase shift pro- 
portional to the local intensity /, which amounts to a 
shift of the scatterer's resonance frequency. For small a, 
Eq. ^ reduces to a x*-^-* -nonlinearity; furthermore, for 
real values of a, the optical theorem is fulfilled, ensuring 
energy conservation. The above set of N complex equa- 
tions ([ij , with the replacement ^ , is formally written as 
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FIG. 1: (Color online) For a given random configuration 
of A'^ = 1000 scatterers with density nA^ — 1, the top plot 
depicts the backward intensity /out as a function of the non- 
linearity strength a. For a > 0.3, several stationary solutions 
coexist. The middle plot shows the smallest module of all 
eigenvalues of the stability matrix M. At each turning point, 
M exhibits a vanishing eigenvalue. The bottom plot depicts 
the smallest real part of all eigenvalues of M. The dashed 
line indicates the stability threshold defined by Re(A) = 0. 
At a ~ 0.1, the solution becomes unstable, and the instabil- 
ity then increases with a. 



F(E, a) = 0, where F and E are real vectors of dimension 
2N. 

Starting from the linear solution (a — 0) for a given 
configuration of the scatterers, a numerical scheme, based 
on a Newton-Krylov method, allows to continuously fol- 
low the solution E(a), which fulfills 



dE 

da 



M 



da 



(3) 



with the Jacobian matrix M = If M has a vanishing 
eigenvalue, this gives rise to a turning point in the curve 
E{a), such that several solutions of Eq. ([T|) coexist for 
the same value of a. 

As an example, we plot in Fig. [T| (top) the intensity 
lout , scattered into the direction opposite to the incident 
wave as a function of a, for a specific configuration of 
N = 1000 scatterers randomly distributed in a sphere 
with density n — N/V such that nX'^ = 1. For a — 0, 
these parameters correspond to mean free path £ = 20/fc 
and optical thickness (measured along the diameter 2R of 
the sphere) 2R/io = 4. In Fig. [T] (top) , we clearly notice 
the existence of several stationary solutions for a > 0.3, 
and in particular for a > 0.6 (see inset). Fig. [TJmiddle) 
shows the smallest module |A|i„in of all eigenvalues of the 
matrix M. At every turning point, the curve |A|mi,i(a) 
touches the horizontal axis, as expected. 

However, even if multiple stationary solutions coexist, 
one will observe speckle instability, i.e. a time dependent 



solution persisting at long time, only if these solutions 
are dynamically unstable. To investigate the dynamical 
behaviour, we use the following relaxation model for the 
time evolution of the dipole d{t): 



d{t)^-^[d{t)~g{I)E) ), 



(4) 



where r/2 corresponds to the inverse relaxation time. E 
and / are now the instantaneous local field and intensity. 
We assume F ^ c/R, such that the field propagation 
can be considered as instantaneous within the radius R 
of the cloud. According to this model, small deviations 
SF,{t) = E(t) -E^* from the stationary solution E^' fulfill 
the following set of equations: 



F 



MSE, 



(5) 



Hence, the stationary solution is stable if all eigenval- 
ues of the matrix M have a positive real part. In Fig. [l] 
(bottom), we plot the smallest real part of all eigenval- 
ues along the stationary solution shown on the top. Just 
above a = 0.1, the solution becomes unstable and the 
instability increases with a. Furthermore, in the region 
of multiple solutions (see around a = 0.5, for instance), 
the different solutions have different instability rates. As 
a closer inspection reveals, the branch between the two 
turning points is the most unstable one, exactly like in 
the standard bistability scheme [TB]. However, whereas 
in the usual scheme the two other branches would corre- 
spond to stable solutions, leading to a bistable behavior, 
they are, in the present case, also unstable. 

In summary. Fig. [TJbottom) reveals that, for the spe- 
cific configuration of scatterers under consideration, no 
stable stationary solution exists for a > 0.1. Hence, we 
expect that, in this regime, the actual dynamics of the 
system does not converge towards a stationary state, but 
remains time-dependent even in the long-time limit. This 
is confirmed by Fig. |2] where we plot, for the same con- 
figuration of scatterers as in Fig. [T] the backscattered 
intensity as a function of time for different values of the 
nonlinearity strength a. One clearly observes a transition 
from the stable stationary solution at a = to a chaotic- 
like behavior at a = 1. For intermediate values, stable 
periodic solutions are obtained, indicating that, for those 
values, the stationary solutions are already unstable. Let 
us stress that the bifurcation from the stable to the un- 
stable regime does not necessarily occur at values of the 
nonlinearity for which many stationary solutions coexist. 
Indeed, since our dynamical system has a very large num- 
ber of degrees of freedom {2N) and, in addition, is not 
a Hamiltonian one, many possible bifurcation scenarios 
are in principle possible, including for example strange 
attractors jlT]. A more detailed analysis would require 
the complete determination of the underlying dynamical 
structure (periodic orbits, limit cycles...), which is be- 
yond the scope of this Letter. 
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FIG. 2: Typical evolution of the intensity in the backward 
direction as a function of time for the same configuration of 
scatterers as in Fig. [l] and for different values of the nonlin- 
earity strength a. One clearly observes a transition from a 
stable stationary solution for a = to a chaotic-like behavior 
for a = 1. For intermediate values, stable periodic solutions 
are obtained, emphasizing that, for those values, the station- 
ary solutions are already unstable. 



The preceding discussion refers to a given single con- 
figuration of scatterers. Although it has been chosen as a 
representative example of the transition to the unstable 
regime, more information can be obtained when perform- 
ing a statistical analysis of the system properties over an 
ensemble of different configurations. Fig. [s] (top) displays 
the cumulative distributions of instability thresholds, de- 
noted Pinst(Q^), for different system parameters. Hence, 
we have plotted, for different numbers N and optical den- 
sities nX^ of scatterers, the proportion of unstable con- 
figurations as a function of the nonlinearity strength a. 
As expected, if we fix the number of scatterers, the in- 
stability occurs earlier (i.e. the distributions are shifted 
towards smaller a) for increasing density (i.e. disorder 
strength). On the other hand, if we fix the density, insta- 
ble configurations are more frequently encountered with 
increasing number of scatterers (i.e. system size). How- 
ever, both variations are not independent: as the figure 
clearly reveals, the statistical distribution depends only 
on a single parameter, namely p = N x nX^. In partic- 
ular, the average threshold, (ainst), depends linearly on 
p~^, see Fig. [s] (bottom, left). 

As we have checked, for large p, the resulting instability 
thresholds are small enough so that first order perturba- 
tion theory with respect to the nonlinearity strength a 
can be applied. Hence, we expect the statistics of insta- 
bility thresholds to be closely related to properties of the 
linear system, in particular the real parts of the eigenval- 
ues of M(q! = 0). The latter define the decay rates F** of 
the linear system [T8j [19] . Since the occurence of a single 
negative eigenvalue is already sufficient for the instability 
to occur, the state with the smallest decay rate (or, equiv- 
alently, longest lifetime) attains a particular importance. 
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FIG. 3; (Color online) On the top, the cumulative distribu- 
tions of instability thresholds are plotted for different numbers 
A'' and densities nA^ of scatterers. A scaling law is revealed, 
according to which the distributions depend only on a single 
parameter p — N x nX^. The bottom left plot shows the av- 
erage threshold to depend linearly on The bottom right 
plot displays the distribution P'^^T") of the smallest linear 
decay rate Fs, which exhibits the same scaling law in terms 
of the parameter p. This strongly emphasizes the crucial role 
played by the interference effects for the speckle instability. 



Therefore, we have plotted in Fig. [sjbottom, right) the 
cumulative distribution of the smallest rate V for dif- 
ferent numbers N of scatterers and densities nX^. Here, 
the solid lines correspond to = 2000 and (from left to 
right) nX^ =0.5, 0.25 and 0.125, the long-dashed lines to 
A^ = 1000 and (from left to right) nX^ =1, 0.5 and 0.25 
and the dotted lines to A^ = 500 and (from left to right) 
nX^ =2, 1. and 0.5. And indeed, we find the distribu- 
tions to be governed by the same parameter p = N x nX^ 
which also determines the statistics of instability thresh- 
olds. From these data, the average value F/ (F*), i.e. the 
longest lifetime, is found, as a crude estimate, to scale 
like p^/^. This time is much longer than the Thouless 
time, i.e the longest time predicted within diffusion ap- 
proximation, which, in F^^ unit, scales like (pnX^)'^^^ 
(i.e. the square of the optical thickness). This differ- 
ence strongly emphasizes that interference effects, which 
enhance the lifetime of a particular single mode with re- 
spect to the diffusive Thouless time, are deeply involved 
in the speckle instability. Due to its long lifetime, such 
a " prelocalized" mode [20] almost behaves like a perfect 
cavity, and is therefore much more severely affected by 
the nonlinearity than short-lived modes and, thus, even- 
tually allows instability to settle. We note that a similar 
argument was also put forward to explain the appearance 
of discrete peaks in the emission spectrum of the coher- 
ent random laser [5T]. Let us note that the instability 
threshold of our nonlinear point scatterers system differ 
from the one predicted in the case of linear scatterers 
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FIG. 4: (Color online) Configuration average of the intensity 
{I{8,t)), scattered at time t = lOOOF"^, in the direction 6, 
where 6 — corresponds to exact backscattering, for 1000 
difi'erent configurations of N = 1000 scatterers with density 
nX^ = 1. Although the coherent backscattering cone at first 
decreases with increasing nonlinearity, the occurence of a dip 
at larger nonlinearities proves that coherent effects survive in 
the time dependent regime a > 0.3, and even in the fully 
chaotic one a > 1. 



ties of the speckle patterns, and observed the transition 
from a regime with a unique, stable and stationary solu- 
tion towards chaotic-like behavior. In addition, we have 
shown that the statistical properties of the instability 
thresholds are determined by the same parameter as the 
smallest decay rates of the linear system, putting thus 
forward the crucial role played by interference effects for 
the speckle instabilities. Furthermore, we have shown 
that, quite surprisingly, the coherent backscattering ef- 
fects are not erased in the time dependent regime, even 
not in the chaotic one. Obviously, it would be very inter- 
esting to perform similar studies in the strongly localized 
regime {k£ < 1). Although our argument based on the 
decay rates of the linear system would then suggest the 
appearance of instabilities for any small amount of non- 
linearity, our results also leave open the possibility to 
observe strong localization even in presence of nonlinear- 
ity, depending on the timescale on which the instabilities 
develop. The authors would like to thank S. Skipetrov, 
C. Miniatura and D. Delande for fruitful discussions. 



in a x^"^-' nonlinear medium for which no instability 
should occur if the nonlinearity is small enough for first 
order perturbation theory to apply. Further research is 
necessary in order to characterize the relevant proper- 
ties of nonlinear disordered systems, leading to different 
scenarios for the development of speckle instabilities. 

Finally, we address the question of coherent transport 
in the unstable regime. For this purpose, we compute 
the coherent backscattering (CBS) cone {I{9,t)), i.e. the 
configuration average of the intensity scattered at a time 
t in the direction 9, where 9 — corresponds to the exact 
backward direction. Numerically, we observe that after, 
typically, few lOOF^^, the average intensity {I{9,t)) be- 
comes time-independent. We plot this stationary value as 
a function of 9 in Fig. [4] where the average is performed 
over 1000 random configurations with the same param- 
eters {N — 1000 and — 1) as for the configuration 
examined in Figs. [T] and [2] Obviously, the net effect of 
the nonlinearity is the reduction of the CBS cone height, 
and the formation of a dip at larger nonlinearities [8j . Ac- 
cording to the theory presented in , this effect can be 
qualitatively explained by a dephasing between reversed 
scattering induced by the nonlinearity. In the present 
context, however, the most important point is that a co- 
herent effect still survives in the time dependent regime 
at a > 0.3, and even in the fully chaotic one at a > 1. 
A possible explanation can be put forward when consid- 
ering the time evolution depicted in Fig |2] the typical 
timescale is much larger than the scatterer response time 
r~^, such that the intensity pattern is almost fixed for a 
photon scattered along short paths inside the medium. 

In summary, using a model of nonlinear point scatter- 
ers, we have analysed the dynamical instability proper- 
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